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We present high-resolution direct numerical simulations of turbulent three-dimensional Rayleigh- 
Benard convection with a focus on the Lagrangian properties of the flow. The volume is a Cartesian 
slab with an aspect ratio of four bounded by free-slip planes at the top and bottom and with pe- 
riodic side walls. The turbulence is inhomogeneous with respect to the vertical direction. This 
manifests in different lateral and vertical two-particle dispersion and in a dependence of the dis- 
persion on the initial tracer position for short and intermediate times. Similar to homogeneous 
isotropic turbulence, the dispersion properties depend in addition on the initial pair separation and 
yield a short-range Richardson-like scaling regime of two-particle dispersion for initial separations 
close to the Kolmogorov dissipation length. The Richardson constant is about half the value of 
homogeneous isotropic turbulence. The multiparticle statistics is very close to the homogeneous 
isotropic case. Clusters of four Lagrangian tracers show a clear trend to form flat, almost coplanar 
objects in the long-time limit and deviate from the Gaussian prediction. Significant efforts have 
been taken to resolve the statistics of the acceleration components up to order four correctly. We 
find that the vertical acceleration is less intermittent than the lateral one. The joint statistics of 
the vertical acceleration with the local convective and conductive heat flux suggests that rising and 
falling thermal plumes are not associated with the largest acceleration magnitudes. It turns out 
also that the Nusselt number which is calculated in the Lagrangian frame converges slowly in time 
to the standard Eulerian one. 

PACS numbers: 47.55.pb, 47.27.te, 47.27.ck 



I. INTRODUCTION 

Turbulent convection is one of the best studied funda- 
mental flows in fluid dynamics research P, Q ■ One reason 
is the large range of examples and applications in nature 
and technology for which a turbulent motion is initiated 
and sustained by heating a fluid from below and cool- 
ing from above. Almost all of these studies have been 
conducted in the Eulerian frame of reference. They were 
primarily focussed to the mechanisms of local [1, |j, d, @| 
and global 0, H, turbulent heat transfer. 

The Lagrangian perspective of turbulence, in which the 
fields are monitored along the trajectories of infinitesimal 
fluid parcels, has recently produced new insights into the 
local topology of fluid parcel tracks, the local strength 
of accelerations, and the statistics of time increments of 
turbulent fields [ll|, ■ The progress is caused on the 
one hand by significant innovations in the experimen- 
tal techniques, such as three-dimensional particle track- 
ing 13, 14, 15] or acoustic methods pji]. On the other 
hand, direct numerical simulations of turbulence become 
now feasible that resolve three-dimensional Lagrangian 
turbulence at moderate and higher Reynolds numbers 
[l3, m, m, HOl ■ Both, experiments and simulations, made 
a deeper understanding of the small-scale intermittency 
and its connection with large accelerations of fluid parcels 
possible. 

Lagrangian investigations in convective turbulence are 
however rare. Several reasons can be given for this cir- 
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cumstance. First, on the experimental side it is de- 
sirable to monitor the temperature along the particle 
tracks beside the velocity components and the acceler- 
ations. Only recently, Gasteuil et al. [2l| constructed 
therefore a smart particle, that monitors velocity, tem- 
perature and orientation while moving through the cell. 
Due to the integrated power supply the particle diam- 
eter remained however larger than the thermal bound- 
ary layer thickness, such that the large-scale bulk motion 
can be monitored only. Second, it is also clear that the 
complexity of direct numerical simulations increases since 
the temperature fleld has to be advected in addition to 
the velocity. Temperature tracking along the tracer po- 
sitions requires additional interpolations. Furthermore, 
one cannot return to simulations in a fully periodic cube, 
the so-called homogeneous Rayleigh-Benard convection 
setup, since the periodicity in the direction of the mean 
temperature gradient causes a self- amplifying fluid mo- 
tion. This was discussed in detail by Calzavarini et al. 
[13, HI]- Third, the turbulence is inhomogeneous -at least 
in the vertical direction as in the following setup- and it 
is thus not clear which of findings from the homogeneous 
isotropic purely hydrodynamic turbulence pertain. For 
example, the height dependence of the statistics has to 
be considered additionally. 

First numerical attempts have been made recently to 
study some aspects of the heat transfer and tracer dis- 
persion in the Lagrangian framework of convective tur- 
bulence [131 • The motivation of the study can be con- 
densed in one question: Which new insight into the na- 
ture of turbulent convection provides the complementary 
Lagrangian view? One result of [13 was to determine 
a mixing zone which is dominated by rising and falling 
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thermal plumes. This is done by combining accelera- 
tion and local convective heat flux statistics. The mixing 
zone starts right above the thermal boundary layer and 
extends several tens of the boundary layer thickness into 
the bulk of the cell. Thermal plumes are fragments of 
the thermal boundary layer that detach in the vicinity 
of the top and bottom isothermal planes. The existence 
of a mixing zone has been suggested in several Eulerian 
studies on the basis of other criteria, e.g. 0,11^ and was 
thus confirmed in the complementary Lagrangian frame 
of reference [l^l • 

The present work extends the previous study [2^ into 
several directions. Beside the local convective, the local 
conductive heat flux is studied along the tracer tracks. 
It requires to monitor temperature gradient components. 
Furthermore, the analysis of the Lagrangian tracer dis- 
persion is extended. In addition to the hydrodynamic 
case , we study the dependence of pair dispersion on 
the initial separation and the initial seeding position. As 
discussed in Refs. [l^,[23,|2l] for the pure hydrodynamic 
case, higher order particle statistics requires to track lit- 
tle clusters of tracers. We provide here an analysis of 
the four-particle-statistics, where the tracers start out of 
groups of tetrahedra of different sidelengths and initial 
vertical positions. 

The outline of the manuscript is as follows. In the next 
section the equations of motion, the numerical scheme, 
the Lagrangian tracer tracking and the turbulent heat 
transfer. In section III, some results of the Eulerian 
statistics of the temperature field are presented. This 
section is followed by sections on the Lagrangian particle 
dispersion, the acceleration statistics and the conductive 
and convective heat flux. We conclude with a short dis- 
cussion of our results and will give a brief outlook to 
possible extensions of the present work. 



II. NUMERICAL MODEL 



a. The temperature field is decomposed into a linear 
profile and fluctuations 9 about the profile 

AT 

T{x,t) = - — {z-H/2) + e{x,t). (4) 

Since T is prescribed and constant at bottom and top 
boundaries z = and z = H, the condition 9 = follows 
there. Here, AT > 0. The dimensionless control para- 
meters are the Prandtl number Pr, the Rayleigh number 
Ra, and the aspect ratio F, 



Pr 
Ra 
F 



agH^AT 
L 

H ' 



(5) 
(6) 
(7) 



The simulation domain is V = LxLxH = [0,F7r] x 
[0,F7r] X [0, tt]. In lateral directions x and y, periodic 
boundary conditions are taken. In the vertical direction 
z, free-slip boundary conditions are used which are given 
by 



and dzUx = dzUy 



0. 



(8) 



The computational grid has a size of Nx x Ny x Nz = 
2048 X 2048 x 513 points. For an aspect ratio F = 4, 
it is thus equidistant in all three space directions with 
a grid spacing Ax. Time-stepping is done by a second- 
order predictor-corrector scheme. The production runs 
are conducted on one rack of the Blue Gene/P system 
which corresponds with 4096 MPI tasks [l^. We use 
volumetric Fast Fourier Transforms based on the pSdfft 
package by D. Pekurovsky [30]. The spectral resolution is 
kmaxTlK = 4.5 where kmax = 2V2TrNx/iSLx). Quantity 
rjK = v^^'^ I (e)^^'* is the Kolmogorov scale with the mean 
energy dissipation rate (e). 



A. Equations of motion and boundary conditions 

The Boussinesq equations, i.e. the Navier-Stokes equa- 
tions for an incompressible flow with an additional buoy- 
ancy term ag9ez and the advection-diffusion equation 
for the temperature field, are solved by a standard pseu- 
dospectral method for the three-dimensional case [29| . 
The equations are given by 



V • M = 



0, 



kV'^9 + Uz 



AT 

IT 



(1) 

(2) 
(3) 



Here, u is the turbulent velocity field, p the (kinematic) 
ure field and 9 the temperature fluctuation field. 



pressure beld and f/ the temperature fluctuation 
The system parameters are: gravity acceleration g, 
matic viscosity i/, thermal diffusivity k, vertical temper- 
ature gradient AT/ and thermal expansion coefficient 



B. Lagrangian particle tracking 

Lagrangian tracer particles follow the streamlines of 
the turbulent velocity field in correspondance with 



X = u{x{t), t) . 



(9) 



For the majority of the analysis, we seeded 3 x 10^ tetra- 
hedra aligned along the outer coordinate axes in the box, 
i.e. Xi = Xq, X2 = Xq + isx, X3 = Xq + ley, and 
X/^ = Xq + Ibz ■ The vector Xq is randomly chosen in the 
box. The tracer ensemble was divided into 6 six groups 
with initial sidelengths of ^ = 1,2,4,8,16 and 32 grid 
spacings Aa; which correspond with 0.5, 1, 2, 4, 8 and 16 
TjK- The two-particle dispersion analysis is consequently 
conducted for the three tracer pairs {xi,X2}, {xijX^}, 
and {xi,X4} of each tetrahedron. 

For the two- and multiparticle statistics, we run in ad- 
dition a simulation with the following initial conditions: 
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FIG. 1: (Color online) Side view of two instantaneous La- 
grangian tracer distributions ai t = 9.5r,, (top) and t = 19t^ 
(bottom). The tracers are colored with respect to the local 
temperature T at their position. Four intervals are taken: 

T G [Ar/2,Ar/3], r e [Ar/3,Ar/6], r e [Ar/6,o] and 

for T e [0, -Ar/6]. The whole ensemble consists of 2 x lO'^ 
tracers. Tracers started in the whole x — y plane at the height 
of the thermal boundary layer thickness which is indicated by 
the solid line very close to the bottom plane in both plots. 




Mean temperature profile 



FIG. 2: (Color online) Mean temperature profile {T(z))A,t of 
the turbulent convection run at i?a = 1.2 X 10* and Pr = 0.7. 
The inset shows the resolution of the thermal boundary layer 
with 7 grid planes. It also indicates the geometric interpre- 
tation of the thermal boundary layer thickness St (see the 
text). The lower right box in the main figure indicates the 
size of the magnification. 



selt number which is given for a plane at fixed height z 

by 



again Xi = Xq, X2 = xq + ie^, x^ = Xq + £ey, and 
X4 = Xq + Iez- The X and y coordinates of the vector 
Xq are again randomly chosen. The vertical coordinate 
corresponds with zq = St/2, 6t, 106t, 20St and H/2. 
Here, we pick £ = riK/2 and 2'rjK- 

The Lagrangian particles are advanced in time simul- 
taneously with the Boussinesq equations. The veloc- 
ity, temperature and temperature gradient components 
at intergrid positions are calculated by trilinear inter- 
polation. The full particle set is written out each 0.45 
Trj. Here, — ^ vj (e) is the Kolmogorov time. Ac- 
celerations along the Lagrangian tracks are calculated 
from three successive integration steps (Ai = 0.006t^) 
and the output interval is the same as for the particle 
positions, velocities, temperature, and temperature gra- 
dient. We thus gather Lagrangian statistics over up to 
4.8 X 10^ tracer particle events. Figure [T] illustrates the 
initial phase of the tracer dispersion. All tracers start 
from a X — y plane close to the bottom wall. 
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where {■)A,t denote averages in planes at z and with re- 
spect to time. The value of Nu{z) is constant and inde- 
pendent of z. The global Nusselt number is then defined 
as 
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H 



Nu = — I Nu{z)dz 
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where {■)v.t is a combined volume and time average. 
The Nusselt number for the present free-slip boundary 
case follows to Nu = 56.36 ± 0.59. Similar to Julien 
et al. (3ll |. we find an enhanced turbulent heat trans- 
port in comparison to no-slip top and bottom plates. For 
Rayleigh numbers between 9.8 x 10^ and 1.2 x 10*, we fit 
the power law Nu = 0.166 x _Ra°-^^^ to the data. 



C. Turbulent heat transfer 

The convective turbulence is studied for one parame- 
ter setting. The Rayleigh number is Ra = 1.2 x 10*, the 
Prandtl number Pr = 0.7 and the aspect ratio F = 4. 
The response of the system is a turbulent heat trans- 
port as quantified by the dimensionless (Eulerian) Nus- 



III. EULERIAN TEMPERATURE STATISTICS 

Figure [5] displays the mean temperature profile as a 
function of height. The total temperature can take values 
between —AT/2 and AT/2 only. As typical for higher 
Rayleigh numbers, the jump of mean profile to zero is 
observed across a thin layer, the thermal boundary layer. 
The inset magnifies the vicinity of the bottom plate. The 



FIG. 3: (Color online) Contour plot of the instantaneous tem- 
perature field T{x,to). 
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FIG. 5; (Color online) Statistics of the temperature fluctua- 
tions w — T' /Trms- Top: Probability density function (PDF) 
of the temperature fluctuations. Data are compared with a 
Gaussian distribution (dashed line). Inset: Skewness and flat- 
ness of the temperature fluctuations w as a function of the ver- 
tical coordinate z/H. The profiles are obtained by averaging 
over lateral planes and a sequence of statistically independent 
snapshots. The Gaussian value for the flatness _F = 3 is indi- 
cated by the dashed line. Bottom: PDFs of the temperature 
fluctuations taken in different planes which are indicated in 
the legend and by vertical solid lines in the inset of the top 
figure. 



FIG. 4: (Color online) Contour plot of the instantaneous ther- 
mal dissipation rate field x(a^,to). Data correspond to those 
in Fig. |3l In order to highlight the small-amplitude dissi- 
pation filaments in the bulk, we plot contours of the decadic 
logarithm of X- 



thickness of the thermal boundary layer is defined as 

H 



6t — 



2Nu 



(12) 



For z = 0, the conductive part of (flOl) contributes to Nu 
only and we can set Nu = ~ndz{T) A,t\z=o/ [n^T / H). 
This leads to 5t = -/^T / {2d,XT) A,t\z=o) and to the geo- 
metric derivation of the thermal boundary layer thick- 
ness (as indicated in the inset of Fig. In con- 
trast to the no-slip case, we have {ux)A.t — {uy)A,t — 
{uz)A,t — . Consequently, no velocity boundary layer 
is present. The Taylor microscale Reynolds number 
Rx = yJh/iM^)) {uD « 143. 

Figure [3] shows an instantaneous snapshot of the to- 



tal temperature field T{x,t). Contour plots in two side- 
planes and close to the top and bottom planes are shown. 
We observe a typical feature of thermal convection - the 
ridge-like maxima which correspond with thermal plumes 
that detach randomly. They form a skeleton which is ad- 
vected by the flow close to the boundaries. The plumes 
coincide with local maxima of the thermal dissipation 
rate field (see Fig. |4] and compare it with Fig. [3]) which 
is defined as 



(13) 



The definition contains the temperature fluctuations 
which are given by 



T'ix,t)^T{x,t)- {Tiz))A,f 



(14) 



The probability density function (PDF) of T' is shown 
in Fig. \5\ We compare the PDF of data taken from 
the whole slab volume with the Gaussian statistics in the 
top figure. Similar to findings for turbulent convection in 



5 



closed cylindrical vessels with solid walls the temperature 
field statistics deviates from Gaussian Q. The analysis 
can be refined. The inset of the top panel shows therefore 
vertical profiles of the plane- and time- averaged flatness, 
F = {T"^)A,t/{T'^)Xt. The flatness differs clearly from 
the Gaussian value of 3 in all parts of the convection cell. 
In addition we plot the profile of the plane- and time- 
aver ae ed skewness S = {T'^) A,t/ {T'^f/^t in the same 
inset. The magntiude of the skewness peaks at about 56t 
which is well inside the plume mixing zone (23 |. Both 
proflles agree also qualitatively with those by Kerr Q 
and by Ref. @] which have been conducted with no-slip 
top and bottom boundaries. In the bottom panel of Fig. 
[5] we show the PDF of the temperature fluctuations in 
four different planes (see the legend). The PDF in the 
midplane comes closest to a Gaussian profile. Our data 
suggest that the free-slip boundary conditions lead to 
smaller deviations from Gaussianity compared to the no- 
slip case. It should also be noted that for strong rotation 
of the cell about the vertical coordinate the temperature 
fluctuations are Gaussian for both, no-slip and free-slip 
boundary conditions, as reported by Julien et al. fsj]. 



IV. LAGRANGIAN PARTICLE DISPERSION 

A. Two-particle dispersion 

The Eulerian framework analysis of turbulent convec- 
tion demonstrated already that the flow is indeed inho- 
mogeneous with respect to the vertical direction. Fur- 
thermore, we recall that the Lagrangian tracer motion is 
constrained between z = and H since both walls can- 
not be penetrated. One motivation to study the disper- 
sion in three-dimensional turbulent convection is there- 
fore to verify if the classical Richardson dispersion law 
[3^ can be also observed for the present case. Recall 
that the Richardson dispersion law follows from a solu- 
tion of a diffusion problem which assumes a homogeneous 
and isotropic turbulent state. It states that, given two 
particle tracks, X2(t) and Xi{t) with Xi = {xi, yi, Zi), the 
distance vector R{t) — X2{t) — Xi{t) will follow 



{R\t))i 



93d 



(15) 



where g^d is a universal constant of 0(1). The symbol 
denotes an average over Lagrangian particle tracks. 
We decompose the relative tracer motion into a lateral 
and vertical contribution in order to separate homoge- 
neous and inhomogencous directions. The distance vec- 
tor can be written as 



(16) 



The lateral two-particle dispersion is given by {R^y{t) — 
R'^y{0)) L where the average is taken over 6 x 10^ particle 
pairs. Here, 

Rxy = [x2{t) - xi{t)]ex + [y2{t) - yi{t)]ey . (17) 
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FIG. 6: (Color online) Vertical and lateral particle pair dis- 
persion as a function of the initial pair separation, (a) Com- 
pensated lateral dispersion for different initial separations as 
indicated in the legend. The dashed lines correspond with 
Eq. (US} for 7?o = ?7if/2 and i?o = r\K and with Eq. ([T9)) 
for i?o = ^r\K and _Ro = (b) Compensated vertical dis- 

persion. Again, Eqs. (|18|l and (|19p are fitted to the data for 
the same initial separations, (c) Same data as in (b) without 
compensation by (e)t^. The vertical dispersion is constrained 
between the planes and levels thus off at larger times. The 
square of the cell height, i?^, and the Kolmogorov time scale, 
r,,, are indicated. All axes are given in decadic logarithm. 
Tracer pairs are seeded initially across the whole volume. 



Similarly, the vertical dispersion is given by {Rl{t) — 
RI{0))l- The dispersion in each space direction would 
contribute with a weight of 1/3 in homogeneous isotropic 
turbulence. In order to compare our pair dispersion re- 
sults with the predictions for isotropic turbulence, we will 
introduce two weight factors, Cxy = 2/3 for the lateral 
motion and Cz = 1/3 for the vertical one. 

Figure [S] displays both dispersion processes with re- 
spect to time for six different initial pair separations as 
explained in section II C. The two-particle dispersion is 
given by a compensated plot in panels (a) and (b) of the 
figure. The graphs are normalized by {e)t^ to capture a 
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Richardson-like scaling as a plateau. The initial ballistic 
behavior at small separations causes then an algebraic 
decay with t^^. Following Sawford et al. fS^I, we fit the 
following two relations to our data at small times 



if i?.m(0) <C rjK and 

{RiM-RimL 
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if Vk ^ Rm{0) ^ L. Here L is the outer scale of tur- 
bulence and C ~ 2 [2^. Index m stands for the lateral 
terms, xy, or the vertical term, z. The agreement with 
(fT8| for the initial Kolmogorov and sub-Kolmogorov sep- 
arations is reasonable. For larger initial separations we 
use (|19p . The larger the initial separation the better 
agree prediction and data. We fitted the two smallest 
and largest initial separations only. None of the initial 
separations is neither much smaller nor much larger than 
the Kolmogorov scale which explains the slight deviations 
of the numerical results from the law s (fT8|) a nd (fT9|) . 

As discussed for example in Refs. [20, [S^], the estab- 
lishment of a Richardson-like regime depends sensitively 
on the initial separation between the tracers. Indeed, for 
one of the six different separations the lateral dispersion 
curve passes through a small plateau with a Richardson 
constant gxy ~ 0.25. This is observed for an initial sepa- 
ration of Rxy{0) = 2riK, (see solid line in Fig. [DJa)). The 
re-translation of the proportionality constant g^y into a 
three-dimensional homogeneous isotropic turbulence case 
is obtained by 



9x^ 



^9xy 



0.375 . 



(20) 



The proportionality constant is smaller than the value 
om ~ 0.5 — 0.6 for homogeneous isotropic turbulence [H, 
IMIlOl. In Ref. [Hi, it was already shown that the PDF 
of the lateral particle pair distance can be fitted to the 
stretched exponential form of Richardson [s^, however 
not to the Gaussian shape as suggested by Batchelor [s^l . 

Figures [6l^b) and (c) display the vertical dispersion. In 
panel (b), we repeat the compensated plot of panel (a) 
and show the fits to A plateau is observed now 

for initial separations between Irjx and 2r]K- The result- 
ing constant is gz ~ 0.05 (see solid line in Fig. (H^b)). 
If one combines the lateral and vertical dispersion, the 
Richardson constant for the turbulent convection follows 
to 



gxy + gz~ 0.3 , 



(21) 



which is less than our earlier estimate of gxy ~ 0.375 
and g^d ~ 0.5 — 0.6. A smaller Richardson constant cor- 
responds with a stronger correlated pair motion. Such 
behavior can be attributed to the presence of rising 
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FIG. 7: (Color online) Vertical and lateral particle pair dis- 
persion as a function of the initial vertical seeding position zq. 
(a) Compensated lateral dispersion for different initial heights 
zo as indicated in the legend (holds for all four figures). The 
initial tracer separation is Rq = '7/f/2. (b) Compensated ver- 
tical dispersion. The initial separation is also Ro ~ tjk /2. 
(c) Compensated lateral dispersion. The initial separation is 
now Ro — 2riK- (d) Compensated vertical dispersion. Again 
Ro ~ 2rjK- All axes are given in decadic logarithm and the 
Kolmogorov time scale is indicated by a dashed line. The 
solid line follows ~ in all plots. 



and falling thermal plumes - a feature that is absent in 
isotropic turbulence. Additionally, it is known that the 
plumes can cluster and form a large scale circulation [3l . 
Figure ^c) demonstrates that the vertical dispersion is 
constrained by the top and bottom planes. The verti- 
cal contribution (i?^(t) - i?^(0))L to the pair dispersion 
levels off. Eventually the lateral dispersion contributes 
solely to the long-time behavior. 

The specifics of the present inhomogeneous flow is that 
not only the initial pair separation, but also the initial 
seeding position is important. This brings us to the sec- 
ond series of particle dispersion studies where tracer pairs 
with fixed distance in different horizontal planes of the 
slab are seeded (see section II C for details) and shorter 
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simulations for about half the duration are rerun. Figure 
[7] summarizes our findings. We picked five intial seeding 
heights: two in the boundary layer, two in the plume 
mixing zone [l^] and the center plane. While panels (a) 
and (b) are for Rq = rix/'i,, panels (c) and (d) are for 
i?o = 277k . The latter is the separation that yielded a 
short Richardson-like range in Fig. Elja). 

Figures El^a) and (c) show that the lateral dispersion 
curves of the tracer subgroups differ in magnitude. The 
local slope is however nearly the same for all subsets. The 
situation is slightly different for the vertical dispersion: 
while the seeding in the center plane causes a gradual 
variation of the local slope of the dispersion curve (see 
Figs. [7]Jb) and (d)), the seeding in the thermal boundary 
layer leads to significant differences after the initial bal- 
listic period. The same result is observed in Fig. [Tfd). 
The reason is that the tracer pairs probe then the detach- 
ment of the boundary layer fragments to full extent. This 
is not the case when starting in the bulk of the cell, ft 
can also be observed that a plateau (which would imply 
Richardson-like scaling) depends sensitively on the intial 
separation and seeding height. 

To summarize this part, Richardson-like dispersion ap- 
pears for a very small range of scales in the present flow. 
Similar to previous studies, we confirm that the initial 
pair dispersion depends sensitively on both, the initial 
pair separation and the initial vertical position of a tracer 
pair in the volume. The qualitative behavior of the tracer 
dispersion in the convection flow is very similar to that in 
homogeneous isotropic turbulence, given the same range 
of Reynolds numbers. The specifics of the turbulence, 
such as the particular driving mechanism or a present 
inhomogeneity, manifests however in the proportional- 
ity constant g^d and causes eventually quantitative de- 
viations from homogeneous isotropic turbulence. Figure 
[Tljc) illustrates this fact very nicely. The scatter of the 
plateaus can be interpreted as a measure of the sensitiv- 
ity. 



B. Multiparticle statistics 




FIG. 8: (Color online) Time evolution of the eigenvalues U 
(see Eq. (|25p for their definition) of one particular four- 
particle cluster. 



mation of the cluster at small times. It is however open, 
what will be observed in the long-time limit. 

The particle tracks Xi{t)^ (i), x^{t)^ and Xi{t) can 
be transformed into the center-of-mass coordinate 



r{t) 



I 4 



(22) 



and the three relative coordinates (which are of interest 
here) 



^[X2{t) ~ Xi{t)] 
^[2x^{t) - X^{t) - X2{t)] 

[ix^it) - xi{t) - X2{t) - xsit)] . (23) 



P2(t) 
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The radius of gyration follows in this frame to Rg = 

\/ p\+ p\+ p\- The shape evolution of the particle clus- 
ter is monitored by the following moment-of-inertia ten- 
sor 



The Lagrangian statistics of higher-order moments re- 
quires to follow more than two Lagrangian tracers si- 
multaneously. In the following, we will focus to the four- 
particle case. The tracers are initially seeded at the edges 
of tetrahedra as discussed in section II C. The distortion 
of such a small particle cluster by the turbulence has 
been studied for the pure hydrodynamic case in Refs. 
[Tgl . [27I [28I . [35} . The original motivation for such analysis 
was to get a deeper geometrical insight into the formation 
of front-like structures in scalar turbulence: in the vicin- 
ity of steep scalar gradients small particle clusters be- 
come co-planar. Furthermore, since the cluster evolution 
probes the whole range of scales of turbulence, one hopes 
to disentangle systematically correlated large-scale ad- 
vection from decorrelated small-scale motion. The pres- 
ence of thermal plumes in convection will alter the defor- 



1=1 



1 1 



(24) 



where a,b = x,y, z is the component index of the vector 
Pi and i — 1,2,3. The real eigenvalues ffi > 32 ^ 53 ^ 
quantify the shape of the particle cluster. Isotropic ob- 
jects correspond with 51 = 52 = 33, cigar-shaped clusters 
with ^ 52 ~ 93 and pancake-shaped clusters with 
.91 ~ 52 S> (?3. Figure [5] shows the time evolution of the 
three eigenvalues for one specific 4-particle cloud. The 
eigenvalues are normalized and given by 



Ik 



9k 



T 



for fc=I,2,3. 



(25) 



Thus < /fc < 1. One can observe, that the eigen- 
value variations become smoother with increasing time. 
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FIG. 9: (Color online) Time evolution of the normalized 
eigenvalues of the moment-of-inertia tensor. The average 
is taken over the whole ensemble of tetrahedra. For times 
t > 10(= 87r^) the data converge to values {I3}l 0.84, 
{I2)l 0.15, and — > 0.01. The tetrahedra are seeded 
in the whole volume initially. 




t 



FIG. 10; (Color online) Time evolution of the normalized 
eigenvalues of the moment-of-inertia tensor. The average is 
taken again over the whole ensemble of 5 x 10* tetrahedra for 
each case. Four tracer particles form initially a tetrahedron 
with a sidelength of t^k/S. Three ensembles started in planes 
with S, IOSt and H/2 as given in the legend. 



et al. [36| for homogeneous isotropic turbulence. In Ref. 
[T^ . tetrahedra were excluded that had two points too 
close or too far of each other such that their values are not 
directly comparable with the present ones. Recent three- 
dimensional particle tracking experiments by Liithi et al. 
[S^l report an (/2)l which is also close to 0.16. 

Since the large scales are probed in the long-term limit 
by the particle clusters, this agreement suggests that con- 
vective and isotropic turbulence on this scale do not dif- 
fer significantly. The relative motion within a cluster is 
insensitive to whether the tracer particles are swept by 
large vortex structures or by large-scale circulation. Con- 
sequently, our findings suggest that the constrained verti- 
cal motion (and thus the inhomogeneity) is not important 
in the diffusive long-time limit of the cluster dynamics. 

One can expect that for times t > lO^r,,, the tracers ad- 
vance independently of each other. The present long-time 
means are compared with the result of a joint Gaussian 
distribution of the relative coordinates, p{pi, p2, ps) ^ 
exp[— (p^ + + ^3)]- This ansatz results to (/i)g = 0.75, 
(^2)0 — 0.22, and (13)0 — 0.03 for the three-dimensional 
case which is obtained by Monte-Carlo simulations [2^. 
The reported mean of (12) l = 0.15 is smaller than the 
Gaussian value. 

Figure [TO] reports the dependence of the shape evo- 
lution from the initial position zq. The effects remain 
small, but systematic. The closer the starting position of 
the tetrahedron to the boundary plane, the faster it con- 
verges into the final quasistatic state. Again the stretch- 
ing and deformation is most efficient when the particle 
cluster passes through the mixing zone right above the 
tnermal boundary layer. 

Figure [TT] displays the PDFs of Ik for five instants t > 
87t,,. The plots highlight two aspects. First, there is still 
a big variety in the amplitudes of individual Ik although 
their means remain almost unchanged. Second, a very 
slow drift in the tails is present which is indicated by the 
arrows in Fig. [Sfa) and (c). 



A convergence of the cluster to an almost coplanar ob- 
ject is observable for larger times, as quantified by the 
small value of I3 . It will turn out now that this example 
displays a typical long-time behavior. 

In Fig. [HI we show the time evolution of the Lagrangian 
ensemble average of the normalized eigenvalues. Data 
for different initial sidelengths of the tetrahedra are com- 
pared. The tetrahedra are seeded across the whole vol- 
ume. Similar to the two-point measure, the initial de- 
formation of the clusters depends sensitively on the side- 
length of the tetrahedron. The smaller the initial side- 
length the stronger the initial stretching of the cluster to 
a cigar-shaped object (for t < 15r,,). After about 87r^ 
all curves collapse and the mean values remain almost 
unchanged. Our data yield = 0.84, (/2)l = 0.15, 

(/s)^ = 0.01. Suprisingly, the obtained mean values are 
very close to the findings of Biferale et al. [l^ and Hackl 



V. ACCELERATION STATISTICS 

The top panel of Fig. [T2| shows the PDFs of the three 
acceleration components. Each component is given in 
units of the corresponding root-mean-square value. As 
expected, the distributions of the two lateral compo- 
nents coincide. Table 1 provides the quantitative de- 
tails of the acceleration statistics and lists for exam- 
ple the skewness S{ak) = (afc)/(o^)'^^^ and the flatness 
F{ak) = (o-t) / i'^k)'^ with k = x,y or z. The numbers for 
the lateral components are almost identical. The vertical 
acceleration component has a smaller flatness which is in 
line with a sparser tail of the corresponding PDF. The 
bottom panel of the same figure provides the statistical 
convergence test of the fourth-order moments where the 
product w'^p{w) is plotted vs. w with w — ak/ak,rms- 
It reflects the fundamental difficulty to gather reliable 
statistics for higher-order moments in Lagrangian turbu- 
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FIG. 11: (Color online) Probability density functions of the 
normalized eigenvalues Ik of the moment-of-inertia tensor (a) 
Ii, (b) I2, and (c) /a. The time instants at which the data 
have been analysed are for t > 87t,, as given in the legend 
in (a). The vertical dashed lines mark the long-time averages 
{IkjL- The arrows in the upper and lower panel indicate that 
there is still a slight drift in the tails although the means in 
Fig. [9] remain nearly unchanged. The tetrahedra are seeded 
again in the whole volume initially. 
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FIG. 12: (Color online) Probability density function (PDF) 
of acceleration components with k — x,y and 2. Upper 
panel: PDF plots. Each component is normalized by its cor- 
responding root-mean-square value. Lower panel: Statistical 
convergence test for the 4th moment of the acceleration. The 
total number of events which as been included for the analysis 
are 4.4 x 10®. 



ma^ nun^ g^^^) ^(^^^ 

1.17 1001 -574 -0.093 63.4 (± 16) 

ay 1.16 559 -540 -0.156 64.4 (± 16) 

1.10 518 -680 -0.119 30.3 (± 11) 



lence. Recall that this analysis is conducted over a set of 
4.4 X 10* events. The area which is occupied by the scat- 
ter of the graphs in the tails of the PDFs determines the 
error bar of the fourth-order moment (and consequently 
of the flatness). We checked that the second and third- 
order moments display almost no scatter (not shown). 
The issue of statistical convergence has been discussed 
for turbulence measurements in a swirling flow [bsI . [39j 
and in numerical simulations of homogeneous isotropic 
turbulence Our values for the flatness F{ak) of the 
lateral flatness are of the same magnitude as those re- 
ported in [38| . 

Figure [T3l refines the statistical analysis of the acceler- 
ation components. Due to the vertical inhomogeneity, we 
report the height dependence of the acceleration statistics 
for one lateral and the vertical component, respectively, 
and plot contours of the joint PDF p{ai, z). The largest 
lateral accelerations and the fattest tails are found close 
to the top and bottom planes. It will turn out in the next 
section that the vorticity is concentrated in cyclones and 
anti-cyclones close to the thermal boundary layer which 
can rationalize the large lateral accelerations. The sup- 



TABLE I: Root-mean-square values, total maxi- 
mum/minimum amplitudes, skewness S{ak) = {a\) / {a\)^^'^ 
and flatness F{ak) ~ (a^)/(a|)^ of the acceleration com- 
ponents. The error bars of F{ak) have been obtained by 
measuring the area of the scatter in the lower panel of Fig. 
1121 Maxima and mininma are given in units of the gravity 
acceleration g. 



port of the PDF decreases monotonically to the center 
plane. We will get back to this point later in the text 
when discussing the role of the vertical vorticity compo- 
nent in connection with plume detachments. In contrast 
to the result for the lateral accelerations, the support of 
the joint PDF p{az, z) shows no significant variation with 
height. It shrinks to zero in the boundary planes (since 
is Uz = 0) and grows rapidly up to about the thermal 
boundary layer thickness. The slight asymmetry of the 
inner contour lines (for the largest probability density 
levels) corresponds with the rising plumes which detach 
from the bottom plane at 2: = and have > and with 
falling plumes a,t z = H for which < 0. Note also that 
the support of all pdfs is the same in the center of the 
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FIG. 13: (Color online) Joint probability density function of 
acceleration components ai and the height z. Left: compo- 
nent ax- Right: component a^. The contours are displayed 
units of the decadic logarithm. The normalization of the ac- 
celeration component in both panels is as in Fig. 1121 
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FIG. 14: (Color online) Time traces of the convective (top) 
and conductive (bottom) heat transfer along two of the 1.2 x 
10® tracers. The Eulerian value of Nu is indicated as a dotted 
line in the top panel. 



cell. This is consistent with the idea that the turbulence 
is close to isotropic far away from the isothermal walls. 

The important result of this section is that there is dif- 
ferently strong intermittency for the vertical and lateral 
accelerations in thermal convection. It is caused by the 
higher level of intermittency in and close to the thermal 
boundary layer. As a consequence, we will have to take 
a closer look at the mechanisms of local heat transfer in 
the vicinity of the top and bottom planes. This is done 
in the next section. 



VI. LAGRANGIAN CONVECTIVE AND 
CONDUCTIVE HEAT FLUX 

A. Lagrangian heat flux 

The local heat flux contributions can be probed in the 
Lagrangian frame of reference. We adopt therefore defi- 
nition PH)) and calculate the local Lagrangian conductive 
and convective flux contributions. They are given by 



Nu{x) = NUconv{x) + NUcond{x) 



H 
kAT 



{u,{x)T{x) ^ Kd,T{x)) , (26) 



along the Lagrangian tracks x{t). Figure [HI shows typi- 
cal time traces of Nu^onv and Nucond along two tracers. 
They display a big variability with respect to time. Even 
negative values are possible for both contributions. As 
expected, the convective term has a significantly larger 
magnitude than the conductive one. The latter is dom- 
inant in the thermal boundary layers where the thermal 
dissipation rate x has the largest magnitude. Figure [15] 
displays the PDFs of the convective and conductive con- 
tributions gathered along the tracks of the whole tracer 
ensemble. In addition, we display the results for the prod- 
ucts UxT and UyT. All quantities are shown in units 
of their root-mean-square values. While the PDFs for 
UxT and UyT are symmetric, those of u^T and ndzT are 
strongly skewed. This reflects the vertical net transfer 
of heat through the volume. The Lagrangian average 
Nul = {Nu{x)) results to 



Nuj 



1 



H 
kAT 



L,t ■ 



(27) 



We have directly verified from the corresponding PDF in 
Fig- [TSTc) that the mean of the Lagrangian conductive 
heat transfer is with 1.04 very close to one. Further- 
more, it is found that Nu^ < Nu. This is in contrast 
to the experimental findings for the smart particle probe 
of Gastcuil et al. [2l| where Nu < Nul- The reason 
for this difference might be due to the finite extension 
of the smart particle that exceeded the thickness of the 
thermal boundary layer. We have compared the result 
for the full record (F) with that of a smaller subset (S) 
which is one third of set (F) and the results differed by 
6.5 % only (see also Fig. [151). ^ time-resolved analy- 
sis shows that NuL{t) = {Nu{x,t))L relaxes slowly to 
the Eulerian value of Nu. The tracers which have been 
seeded randomly or at particular heights at the beginning 
have to pass a kind of "thermalization" process. 

Our result sheds interesting light on the joint velocity- 
temperature sampling properties of the Lagrangian trac- 
ers. First, it is known that for the present geometry a 
large-scale circulations are present 0, |4l| . Tracers will 
preferentially follow the circulation motion in the convec- 
tion layer. The large-scale circulation can carry a frac- 
tion of the total heat transport only, as has been analyzed 
recently with a Proper Orthogonal Decomposition [43 |. 
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FIG. 15: (Color online) Probability density functions (PDF) 
of the conductive and conductive heat flux contributions along 
the Lagrangian tracer tracks, (a) Lateral convective contri- 
butions UxT and UyT. (b) Vertical convective contribution 
UzT. (c) Conductive contribution ~dT/dz. All quantities 
are normalized by their root-mean-square values. The anal- 
ysis is conducted over two different data sets: the full record 
(F) and a third of it (S). 



Second, we observe high-amplitude events of the vertical 
vorticity lJz — dxUy — dyUx close to the maxima of the 
thermal dissipation rate x- This is shown in Figs. [H] 
and [T7] where contour plots of slice snapshots of x and 
Uz at a height z — 6t are compared. The local maxima 
in the dissipation rate plot (fTB|) reproduce the skeleton 
of plume sheets. In their vicinity, we observe cyclones 
and anti-cyclones that are generated in connection with 



FIG. 17: (Color online) Contour plot of the vertical vorticity 
component ujz{x,y,z = St, to). Data correspond to those in 
Fig. [13 



the detachment of plume fragments. This observation is 
in line with results in Q for the non-rotating and with 
[Sl'l for the rotating case. Exactly these cyclones and 
anti-cyclones cause the large lateral positive and nega- 
tive accelerations as seen in the PDFs in Fig. [121 Fig. [18] 
plots vertical profiles of means of both quantities. While 
the root-mean-square of the vertical vorticity component 
varies weakly, the thermal dissipation rate is strongly 
peaked in and close to the boundary layer. It is also 
known from studies in homogeneous isotropic turbulence 
[ll| that the Lagrangian tracers are not very frequently 
trapped in the core of such vortex structures. 
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FIG. 18: (Color online) Vertical profiles of the averaged ther- 
mal dissipation rate {x{z)}A,t (solid line) and the root-mean- 
square of vertical vorticity component {u!^(z))A,t (dashed 
line). The latter is divided by a factor of 10. 




FIG. 19: (Color online) Joint probability density function 
of the lateral acceleration (upper panel) and vertical (lower 
panel) acceleration and heat flux contributions. Upper panel: 
w = UzT / {uzT)rms. Lower panel: w = u^T j (u^T)rm3- The 
contour spacing is in decadic logarithm. 



B. Joint statistics of Lagrangian heat flux and 
acceleration 



The thermal plumes detach permantly from the ther- 
mal boundary layer and can be identified as regions in 
which the product w^T' > [111, H^. Here we extend 
this analysis and study the correlations between the ver- 



tical velocity component and the (total) temperature in 
relation to the acceleration. Figure [TO] displays the joint 
statistics of the vertical and lateral accelerations and the 
products of velocity and temperature fluctuations, u^T. 
In order to highlight the statistical correlation between 
the two variables of the joint PDF we divide the joint 
PDF by the two single quantity PDFs, 



n(ai,UzT) 



p{ai,UzT) 
p{ai)p{uzT) 



(28) 



The top panel of Fig. [11] shows the joint statistics for 
ax and u^T. A pronunced maximum at larger accel- 
erations and values of UzT is found. They can be re- 
lated to coherent structures, such as vorticity tubes in 
the bulk of the slab or the cyclones/anti-cyclones in the 
boundary layer. The correlation between vertical ac- 
celeration and the product UzT is weaker. The local 
amplitudes of the joint PDF found at the outer bound- 
aries of the support, at moderate acceleration amplitudes 
and larger values of u^T. In Ref. [l^l, we reported the 
same behavior for p{az, UzT') and identified rising plumes 
{az > 0, UzT' > 0) and faUing plumes (a^ < 0, UzT' > 0). 
Recirculations around rising and falling plumes have to 
form due to the incompressibility of the fluid. They were 
related to maxima in the halfplane UzT' < 0. These are 
only some of the possible scenarios which can be assigned 
to strong correlations in the joint statistics. The firm 
conclusion which we can draw from the present analysis 
and the one in [24i | is that plumes (and therefore the ver- 
tical convective flux events) are not connected with the 
largest vertical accelerations. The detachment of plumes 
is a more gradual process. 

We repeated the joint statistical analysis for the con- 
ductive part, —KdzT(x). Results are summarized in Fig. 
[20] The contour plots display now 



n(a„-9,T) 



p{az, -dzT) 
p{az)p{~dzT) 



(29) 



We keep in mind from ([TII)) and (pS)) that upward con- 
ductive heat transport events have a negative sign, i.e. 
~\dzT\. While the lower panel of Fig. [20] includes 
the Lagrangian data of the whole volume, the upper 
panel of the same figure excludes events in the thermal 
boundary layer up to the beginning of the plume mix- 
ing zones. This zone was identified and studied in 
and starts for the present parameter setting at a height 
H/16. The extended tail in the lower panel can thus 
be related to largest gradients (and thus largest thermal 
dissipation rate amplitudes) in and above the thermal 
boundary layer. The asymmetry to negative for the 
largest gradients can be interpreted as tracer decelera- 
tions which are present when temperature gradients are 
formed. Negative accelerations (or decelerations) seem 
to be frequently related to stagnation-point flow topolo- 
gies, those flows which can steepen the temperature field 
to large gradients. 
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FIG. 20: (Color online) Joint probability density function 
of the vertical acceleration component a^/a^^rms and w = 
—dzT/{dzT)rms- Top panel: Tracer positions in the bulk 
have been considered only. Bottom panel: Tracer positions 
in the whole volume are included. The contour spacing is in 
decadic logarithm. 



be obtained for more complex flows than isotropic homo- 
geneous turbulence. The proportionality constants are 
however different which can be attributed to qualitaively 
different turbulence structures such as thermal plumes 
in convection. They affect the vertical dispersion more 
significantly than the lateral dispersion. 

The inhomogcneity of the convective turbulence mani- 
fests in less intermittent statistics of the vertical acceler- 
ation component compared to the lateral ones. Thermal 
plumes are not coupled with the strongest accelerations. 

We find that the Lagrangian Nusselt number converges 
slowly to the Eulerian value. A closer inspection of this 
point could be done in several steps: first, to disentan- 
gle the large-scale circulation from the turbulent back- 
ground, as done in |42| , and to combine such a study with 
a Lagrangian analysis. Second, the present study indi- 
cates also that the Nusselt number relaxes faster with 
a growing number of tracers, in other words the sam- 
pling of joint velocity-temperature statistics improves. 
Our observation is also related to recent experimental 
and numerical studies at very large Rayleigh numbers 

|43| in which the existence and growth of a so-called 
superconducting core is discussed, which is in line with 
a decreasing importance of the large flow circulation for 
growing Rayleigh numbers. A decreasing importance of 
a coherent large-scale circulation might lead again to a 
faster convergence Nul Nu. More studies of the La- 
grangian frame of high-Rayleigh number turbulence are 
thus necessary. 



VII. SUMMARY AND DISCUSSION 

The focus of the present work was on Lagrangian as- 
pects of turbulent convection. The results can be sum- 
marized as follows. The study of pair and multiparticle 
dispersion yields qualitatively similar results compared 
to the homogeneous isotropic case. Although the scaling 
behavior of second order moments is sensitive to initial 
separations and seeding heights, initial separations close 
to the Kolmogorov length result in a short Richardson- 
like scaling range. Interestingly, we reproduce the same 
long-time limits for the particle cluster shapes as in 
isotropic turbulence despite turbulent convection is in- 
homogeneous. This limit deviates from the Gaussian 
value. Our results suggest that the dispersion laws can 
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